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ABSTRACT 

Context. The Sun is a magnetic star whose cyclic activity is thought to be linked to internal dynamo mechanisms. A combination 
of numerical modelling with various levels of complexity is an efficient and accurate tool to investigate such intricate dynamical 
processes. 

Aims. We investigate the role of the magnetic buoyancy process in 2D Babcock-Leighton dynamo models, by modelling more accu- 
rately the surface source term for poloidal field. 

Methods. To do so, we reintroduce in mean-field models the results of full 3D MHD calculations of the non-linear evolution of a 
rising flux tube in a convective shell. More specifically, the Babcock-Leighton source term is modified to take into account the delay 
introduced by the rise time of the toroidal structures from the base of the convection zone to the solar surface. 

Results. We find that the time delays introduced in the equations produce large temporal modulation of the cycle amplitude even 
when strong and thus rapidly rising flux tubes are considered. Aperiodic modulations of the solar cycle appear after a sequence of 
period doubling bifurcations typical of non-linear systems. The strong effects introduced even by small delays is found to be due to 
the dependence of the delays on the magnetic field strength at the base of the convection zone, the modulation being much less when 
time delays remain constant. We do not find any significant influence on the cycle period except when the delays are made artificially 
strong. 

Conclusions. A possible new origin of the solar cycle variability is here revealed. This modulated activity and the resulting butterfly 
diagram are then more compatible with observations than what the standard Babcock-Leighton model produces. 

Key words. Magnetic fields - Sun: dynamo - Sun: activity - Sun: interior - Methods: numerical 



1. Introduction 

Our Sun is a prime example of a very turbulent and magneti- 
cally active star. Its robust 22-yr activity cycle originates in the 
periodic polarity reversal of the Sun's internal large-scale mag- 
netic field. As a result, sunspots emerge at the solar surface with 
statistically well-defined dynamical and morphological charac- 
teristics. These sunspots are thought to be the surface manifes- 
tations of strong toroidal magnetic structures created at the base 
of the convection zone and which rise through the plasma un- 
der the effect of magnetic buoyancy (IParker 19551 1. This toroidal 
field undergoes cyclic variations along with the poloidal field, 
which flips polarity at sunspot cycle maxima. Even if a robust 
regular activity is easily exhibited in the Sun, a significant mod- 
ulation of both the amplitude and the frequency of the cycle has 
been observed. In particular, periods of strongly reduced activity 
have been revealed, the most famous being the Maunder mini- 
mum between 1650 and 1700, during when no sunspots were to 
be seen at the solar surface ( Eddy 1976) l. Today, understanding 
the origins of such a variability has become crucial. Not only is 
the impact on satellites, astronauts, power grids or radio com- 
munications significant but the climatologists community now 
tends to address more and more the question of the influence of 
solar variability on the Earth climate (e.g. IBard & Frank 20061 
[Lean 2005 ). 

The classical explanation for the cyclic activity of the large- 
scale magnetic field is that a dynamo process acts in the solar 
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interior to regenerate the three components of the magnetic field 
and sustain them against ohmic dissipation. The inductive ac- 
tion of the complex fluid motions would thus be responsible 
for the vigorous regeneration of magnetic fields and for its non- 
linear evolution in the solar interior (see lCharbonneau 2005l and 
IMiesch 2005 1 for recent reviews on the subject). 

Understanding how these complex physical processes op- 
erating in the solar turbulent plasma non-linearly interact is 
very challenging. One successful and powerful approach is to 
rely on multi-dimensional magnetohydrodynamic (MHD) simu- 
lations. In this context, two types of numerical experiments have 
been performed since the 70's: kinematic mean-field axisym- 
metric dynamo models which solve only the mean induction 
equation ( Steenbeck & Krause 19691 IRoberts 19721 IStix 19761 
IKrause & Raedler 1980!) and fuU 3D global models which ex- 
plicitly solve the full set of MHD equations dGilman 1983! 
Glatzmai er 19851 ICattaneo 19991 IBrun et al. 2004l l. Those two 
approaches are complementary and needed since 2D mean-field 
models are limited by the fact that they rely on simplified de- 
scriptions of complex physical processes such as turbulence and 
since the cost of 3D models make it difficult, as of today, to pro- 
vide any reliable predictions concerning the large-scale magnetic 
cycle. 

As far as the first kind of numerical simulations is con- 
cerned, a particular model has been favoured by a part of 
the community, namely the Babcock-Leighton flux-transport 
model first proposed by jBabcock (1961)| and Leighton (1969) 
This model has proved to be very efficient at reproduc- 
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ing several properties of the solar cycle as the equator- 
ward branch of activity or the phase relationship between 
toroidal and poloidal fields. It has thus been extensively 



used since the 90's (Wang & Sheeley 1991 , Durne y 1995 
Dikpati & Charbonneau 1999,|Nandy & Choudhuri 2001| l, even 



to make tentative predictions of the next solar cycle 
( Dikpati & Oilman 2006^ . In this formulation, the poloidal field 
owes its origin to the tilt of active regions emerging at the Sun's 
surface at various latitudes during the solar cycle. It thus relies 
on a diff'erent mechanism than "classical" aO. dynamo models 
for which the poloidal field is generated by turbulent helical mo- 
tions twisting and shearing toroidal field lines inside the convec- 
tion zone. The emergence of tilted active regions at the photo- 
sphere is the result of the complex non-linear evolution of strong 
toroidal structures from the base of the convection zone (CZ), 
where they are created through the Q-eff'ect. It is thus natural to 
rely on 3D MHD simulations of rising toroidal structures in the 
convection zone to gain some insight on the best way to model 
the Babcock-Leighton source term. 

Many models carried out since the 80's relied on the as- 
sumption that toroidal flux is organised in the form of dis- 
crete flux tubes which will rise cohesively from the base of the 
CZ up to the solar surface. These models enabled to demon- 
strate that the initial strength of magnetic field was an impor- 
tant parameter in the evolution of the tube and that the active 
regions tilts could be explained by the action of the Coriolis 
force on the magnetic structure (ID'Silva & Choudhuri 19931 
IFan et al. 19941 |CaUgari et al. 1995) . More sophisticated multi- 
dimensional models were then developed and extended to the 
upper part of the CZ and the transition to the solar atmo- 



sphere (pagara 2004|IArchontis et al. 2005l|Cheung et al. 2007 
[Martinez et al. 2008l l. Computations were then performed to 
study the influence of convective turbulent flows on the dynam- 
ical evolution of flux ropes inside the CZ (IDorchet al. 20011 
[Cline 2003" Fan et al. 20031 Jouve & Brun 2009). In particular in 
[Jouve & Brun (20(39)1 such computations were made in a global 
spherical geometry and thus allowed to assess the combined role 
of hoop stresses, Coriolis force, convective plumes, turbulence, 
mean flows and sphericity on the tube evolution and on the sub- 
sequent emerging regions, along with the usual parameters such 
as field strength, twist of the field lines or magnetic diffusion. 
How can the results of 3D simulations of this presumed ma- 
jor step of the dynamo loop be reintroduced in 2D mean-field 
models? This is the question we are addressing in this work. A 
particular point we retain from these simulations is that the rise 
velocity and thus the rise time of magnetic structures is strongly 
dependent on the initial field strength at the base of the convec- 
tion zone. Indeed, magnetic buoyancy competes with convective 
flows to control the trajectory and the speed of the flux tubes 
while they rise. Since the Babcock-Leighton source term relies 
on this process to regenerate poloidal fields, it may be worth try- 
ing to introduce magnetic-field dependent time delays caused by 
rising toroidal fields in those models and analyze the influence 
on the solar cycle. 

Indeed, time delays built up in solar dynamo models have 
been shown to cause long-term modulation of the dynamo 
cycle and under some circumstances to lead to a chaotic 
behaviour (Yoshimu ra 19781 ICharbonneau et al. 20051 or 
IWibnot-Smith et al. 2005T l. In the framework of mean-field 
dynamo theory, several possibilities have been studied to 
explain the variability of the solar cycle. They mainly fell 
into two categories: stochastic forcing or dynamical nonlin- 
earities. Indeed, as we mentioned above, the solar convection 
zone is highly turbulent and it would be surprising if the 



dynamo processes acting inside this turbulent plasma were 
nicely regular The influence of stochastic fluctuations in the 
mean-field dynamo coefficients has been studied in various 
models (e.g. Hoyng 1988| Ossendrijver& Hoyng 1996] 



IWeiss & Tobias 200 01 



t Charbonneau & Dik pati 2000f. 



Moreover, the dynamical feedback of the strong dynamo- 
generated magnetic fields is likely to be significant enough 
to produce non-linear effects on the activity cycle. A 
number of models have introduced these non-linear ef- 
fects (l^octor 19771 iTobias 19971 IMoss & Brooke 20001 
|Bushby 2006| |Rempel 2006| l and have resulted in the produc- 
tion of grand minima-like periods or other strong modulation 
of the cyclic activity. However, time delays have hardly been 
considered and if they were, they were mainly due to the 
advection time by meridional flow. Indeed, the time-scale of 
the buoyant rise of flux tubes was considered to be so small 
compared to the cycle period (and to the meridional flow 
turnover time) that this particular step was assumed to be 
instantaneous. However, we would like to address the question 
of the influence of magnetic field dependent delays on the cycle 
produced by Babcock-Leighton models and especially on its 
potential modulation. We propose to do so in the present paper. 

This article is organized as follows: Section |2] summarizes 
the results of 3D calculations which will be used in our mod- 
ified 2D mean-field Babcock-Leighton model. The formulation 
of this new model is then presented, with a particular focus on 
the modification of the surface source term. Results of 2D mod- 
els are shown in Sect. |4] and analyzed in the following two sec- 
tions. Sect.|5]and|6]present and study the behaviour of a reduced 
set of ordinary differential equations designed to gain some in- 
sight on the results of the 2D model. We discuss the results and 
conclude in Sect. Q 



2. Summary of the 3D calculations 

In this section, we briefly summarize the results obtained by 
|Jouve & Brun (2009)| They use 3D simulations in a rotating 
spherical convective shell to compute the non-linear evolution 
of a magnetic flux tube inside the turbulent convection zone. 



2. 1 . Model and background 



The computations presented in |Jouve & Brun (2009)| use the 
ASH code which solves the anelastic equations of magneto- 
hydrodynamics in three dimensions and in spherical geome- 
try ( Clune et al. T999llMiesch et al. 2000IIBrun et al. 2004J . The 
first step of these calculations is to trigger the convective in- 
stability in the rotating sphere and let the model evolve until 
it reaches a statistically steady state. A twisted magnetic flux 
tube is then introduced at the base of the modelled convection 
zone, in entropy and pressure equilibrium. The lack of mechan- 
ical equilibrium thus causes the tube to be initially buoyant and 
rise through the convection zone up the the top of the computa- 
tional domain. In this work, the interactions between the mag- 
netic tube-like structure and convective motions were investi- 
gated for different initial field strength and these results were 
compared to isentropic cases (where the convective instability 
is not triggered). We do not go into the details of the simula- 
tions and the results. We just wish to summarise some of the 
findings of these 3D calculations which will be useful for the 
present work. 



Jouve, Proctor & Lesur: Time delays and BL dynamo models 



2.2. The rise velocity 



As shown in |Jouve & Brun (2009)| while it rises, the flux tube 
creates its own local velocity field which may strongly disturb 
the background velocity field, especially when the initial mag- 
netic field intensity is strong compared to that of the equiparti- 
tion field Bgq (defined as the field which energy is equal to the 
kinetic energy of the strongest downdrafts). This explains why 
the tube is more or less influenced by the convective motions as 
it evolves in the CZ. Indeed, if the initial magnetic field on the 
axis of the rope is 6 x 10^ G (coiTesponding to 10 B^g), the tube 
creates a velocity through the action of the Lorentz force which 
dominates against the background velocity. Since a strong up- 
flow is thus created near the tube axis, the rising mechanism is 
very efficient and the tube reaches the top of the CZ in only 3 
to 4 days. In the weaker case (where the initial field strength is 
3 X 10^ G, corresponding to 5 B^g) , the velocity field created by 
the magnetic structure is comparable to the background velocity 
field and the latter is thus able to influence the behaviour of the 
flux rope as it rises, the rise time is in this case of about 12 days. 
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Fig. 1. Rise velocity of tubes introduced at 5 B^g (upper panel) 
and 10 Beg (lower panel) measured as a function of time during 
the evolution of the tubes from the base of the CZ up to about 
Q.9Re. 

Figure [1] shows more precisely the evolution of the rise ve- 
locity of tubes introduced at the two different field strengths 
discussed above. On the upper panel, the rise velocity for the 
5 Bpcy-field is shown and the lower panel shows the evolution 
for the stronger field (10 B^g). We first note the different time 
scales for the tubes to reach the surface, and this is the particular 
feature which will be introduced in 2D mean-field models. This 



figure also shows that the acceleration phase of these two tubes 
are similar, but after the maximal velocity has been reached, the 
behaviour of the two tubes differs. In particular, the weak tube 
gets strongly influenced by the convective motions and magnetic 
diffusion. As a consequence, its rise velocity is significantly re- 
duced as it proceeds through the convection zone. On the con- 
trary, after the strong tube has reached its maximal velocity, it 
keeps making its way to the surface at the same speed and the 
strong decrease shown on the figure is eventually strongly linked 
to the upper boundary condition (the boundary conditions for the 
velocity are stress-free and impenetrable in this case). 

3. Reintroduction in 2D models 

The idea of this work is to design a new 2D mean-field flux 
transport Babcock-Leighton model which takes into account the 
findings of 3D calculations concerning the rise of strong toroidal 
structures from the base of the convection zone up to the surface. 
In this section, we thus present the usual mean-field equations 
and the basic ingredients used in the model, with a particular 
focus on the new Babcock-Leighton surface source term. 

3.1. The mean field Babcock-Leighton model 

To model the solar global dynamo, we use the hydromagnetic in- 
duction equation, governing the evolution of the magnetic field 
B in response to advection by a flow field V and resistive dissi- 
pation. 



dt 



V X (V X B) - V X (t/V X B) 



(1) 



where 77 is the magnetic diffusivity. 

As we are working in the framework of mean-field theory, 
we express both magnetic and velocity fields as a sum of large- 
scale (that will correspond to mean field) and small-scale (asso- 
ciated with fluid turbulence) contributions. Averaging over some 
suitably chosen intermediate scale makes it possible to write 
two distinct induction equations for the mean and the fluctuat- 
ing parts of the magnetic field. A closure relation is then used to 
express the electromotive force in terms of mean magnetic field, 
leading to a simplified mean-field equation. In this work we will 
replace the emf by a surface Babcock-Leighton term (Babcock 
1961; Leighton 1969; Wang et al. 1991; Dikpati & Chai-bonneau 
1999; Jouve & Brun 2007a) as described in details below. 

We work here in spherical geometry and under the assump- 
tion of axisymmetry. Reintroducing a poloidal/toroidal decom- 
position of the magnetic and velocity fields in the mean induc- 
tion equation, we get two coupled partial differential equations, 
one involving the poloidal potential A^ and the other concerning 
the toroidal field B^,. 
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where cr = r sin 6, tj, is the turbulent magnetic diffusivity (dif- 
fusivity in the convective zone), \p the flow in the meridional 
plane (i.e. the meridional circulation), Q the diff'erential rota- 
tion, S{r,9,B^) the Babcock-Leighton source term for poloidal 
field. In order to write these equations in a dimensionless form, 
we choose as length scale the solar radius (^o) and as time 
scale the diffusion time (Rq/t],) based on the envelope diffusiv- 
ity (77,). This leads to the appearance of three control parameters 
Cn = HqRI/t],, Cs = sqRq/t], and R^ = vq/Jq/t/, where Qq, sq, vq 
are respectively the rotation rate and the typical amplitude of the 
surface source term and of the meridional flow. 

Equations |2] and [3] are solved in an annular meridional cut 
with the colatitude e [0, n] and the radius (in dimensionless 
units) r e [0.6, 1] i.e from slightly below the tachocline (r = 0.7) 
up to the solar surface, using the STELEM code. This code has 
been thoroughly tested and validated thanks to an international 
mean field dynamo benchmark involving 8 different codes 
JJouve et al. 2008a . At = and 6 - n boundaries, both A^ 
and B^ are set to 0. Both A^ and B^ are set to at r - 0.6. 
At the upper boundary, we smoothly match our solution to an 
external potential field, i.e. we have vacuum for r > 1. As initial 
conditions we are setting a confined dipolar field configuration, 
i.e the poloidal field is set to sin OjP' in the convective zone and 
to below the tachocline whereas the toroidal field is set to 
everywhere. 



3.2. The standard physical Ingredients and the modified 
surface source term 

We now need to prescribe the amplitude and profile of the vari- 
ous ingredients acting in equations |2] and |3] namely the rotation, 
the diffusivity, the meridional flow and the Babcock-Leighton 
source term. 

The rotation profile captures some realistic aspects of the 
Sun's angular velocity, deduced from helioseismic inversions 
(Thompson et al. 2003), assuming a solid rotation below 0.66 
and a differential rotation above the interface. 



Q(r,0) = Q, + i[l+erf(2:!— ^)] 
2 di 



(D.Eq + cLi cos^ Q + a^ cos'' Q - Qc) 



(4) 



bottom boundary r — 0.6 and thus penetrating a little below the 
tachocline. We take a stream function 
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which gives, through the relation Up 
components of the meridional flow 
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(7) 
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rb) 



(8) 



with rb = 0.6. 

In Babcock-Leighton dynamo models, the poloidal field 
owes its origin to the tilt of magnetic loops emerging at the so- 
lar surface. Since these emerging loops are thought to rise from 
the base of the convection zone through magnetic buoyancy, we 
see that we can directly relate the way we model the Babcock- 
Leighton (BL) source term and the results of 3D calculations 
shown in Sect|2l 

In the standard model, the source term is confined in a thin 
layer at the surface and is made to be antisymmetric with re- 
spect to the equator, due to the sign of the Coriolis force which 
changes from one hemisphere to the other These features are 
retained in our modified version. 

However, the standard term is proportional to the toroidal 
field at the base of the convection zone at the same time, imply- 
ing an instantaneous rise of the flux tubes from the base to the 
surface where they create tilted active regions. The 3D calcula- 
tions showed that the rise velocity and thus the rise time depend 
on the field strength at the base of the CZ, we thus introduce 
in our modified version of the source term, a magnetic field- 
dependent time delay in the toroidal field at the base of the con- 
vection zone. We thus take into account the time delay between 
the formation of strong toroidal structures at the base of the con- 
vection zone and the surface regeneration of poloidal field. 

The modified expression of the source term is thus 



with Q.Eq = 1, Oc = 0.93944, r^ = 0.7, dx = 0.05, a^ = 
-0.136076 and 04 = -0.145713 and where erf is the error 
function. With this profile, the radial shear is maximal at the 
tachocline. 

We assume that the diffusivity in the envelope ;; is dominated 
by its turbulent contribution whereas in the stable interior r\^ <K 
7]^. We smoothly match the two different constant values with 
an error function which enables us to quickly and continuously 
transit from 77c to r]x i.e. 



^{f) 



'/c + - (//t 



nz) 



1 H-erf 



(^) 



(5) 



with 77c = lO'^cm^s-' and d = 0.03. 

In Babcock-Leighton flux- transport dynamo models, merid- 
ional circulation is used to link the two sources of the magnetic 
field namely the base of the CZ and the solar surface. For all 
the models presented in this paper we use a large single cell per 
hemisphere, directed poleward at the surface, vanishing at the 



5(r,0,B^) 



l+erf(l^) 



l_erf(l— i) 
d2 



^ ^ g^(r„g,f-TB) 
Bo 



X cos Q sin Q B^ifc, 6, t - tb) 



(9) 



where ra = 0.95, c/a = 0.01, Bo = lO"* G, with the time delay 
Tb proportional to the inverse of the magnetic energy at the base 
of the convection zone, namely 

TB(6,t)^To/B^(r,,6,tf 

In our 3D simulations, the approximate rise time for a 6 x 
10^ G field is indeed about 4 times that of a 3 x 10^ G field. 
This takes into account the more significant effects of convective 
downdrafts and Ohmic diffusion in the "weak B" case. Note that 
the time delay is then dependent both on space and time. 



Jouve, Proctor & Lesur: Time delays and BL dynamo models 



Finally, we take into account the fact that very strong toroidal 
structures (more than 10^ G) are not influenced by the Coriolis 
force enough to gain a significant tilt as they reach the surface. 
To do so, we introduce a quenching term in the surface source 
which will make the regeneration of the poloidal field less ef- 
fective when the toroidal field at the base of the CZ is strong 
enough. Inversely, when the flux tubes are too weak, they will 
not be able to reach the surface at all, and will not take part into 
the regeneration term for the poloidal field. We thus prevent the 
weakest flux tubes (or equivalently the most delayed ones) to 
take part in the regeneration of poloidal fields at the surface. To 
do so, the source term is set to zero when the delays are above a 
certain threshold value corresponding to a rise time of approxi- 
mately half a solar cycle. 

4. Results 

4.1. The Standard solution 

We first quickly present the solutions of our dynamo equations 
for the standard case which has been extensively discussed and 
used in the community to model the solar cycle (e.g. Dikpati & 
Charbonneau, 1999). 

Figure |2] shows the typical solution of a Babcock-Leighton 
standard dynamo model with a realistic choice of param- 
eters. This is the same case as the reference case of 
[Jouve & Brun (2007a)[ In particular, the maximum latitudinal 
velocity at the base of the convection zone (which is the sig- 
nificant velocity for the advection of strong toroidal fields close 
to the tachocline) is set to about 1 m.s"'. 

The butterfly diagram resulting from this particular model is 
in good agreement with the observed migration of sunspots at 
the surface. In particular, a strong equatorward branch stretch- 
ing from the mid-latitudes to the equator is clearly visible on 
the upper panel of Figure |2] A strong poleward branch also 
appears in the evolution of the toroidal field at the interface, 
which is quite typical of Babcock-Leighton models. Indeed, it 
is the result of the combined action of low diffusivity at the 
base of the convection zone and strong radial shear linked to 
the differential rotation profile. It has been shown that both 
higher diffusivities ( [Dikpati & Charbonneau 1999[ ) and penetra- 
tion of the meridional flow below the core-envelope interface 
( [Nandy & Choudhuri 2001 1 could prevent such a branch to ap- 
pear. However, smce our goal is more to study the effects of time 
delays in BL models than to build a calibrated solar-like solu- 
tion, we adopt this solution as our reference case. 

Nevertheless, it has here to be pointed out that making the 
link between the evolution of toroidal field at the base of the con- 
vection zone and the sunspots location at the solar surface dur- 
ing a magnetic cycle is not straightforward. Indeed, to be able 
to make this link, we have to assume that the rise of toroidal 
structures from the base of the convection zone to the surface 
is radial. However, as [Choudhuri & Oilman (1987)[ first demon- 
strated using the thin flux tube approximation and as [Fan (2008)[ 
confirmed with 3D MHD simulations, initially weak magnetic 
flux tubes are deviated from the radial trajectory and tend to fol- 
low a path which is parallel to the rotation axis while they make 
their way to the solar surface. In the simulations presented here, 
we do not take into account this drift in latitude of weak struc- 
tures but intend to do so in a future work. We thus assume the 
evolution of toroidal field at the base of the convection zone to be 
a good proxy for the sunspot migration at the surface but keep in 
mind that the butterfly diagram may take a different form if the 
latitudinal drift is taken into account. 
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Fig. 2. Time-latitude cut of the toroidal field at the base of the 
convection zone (representing the butterfly diagram) and tem- 
poral evolution of the toroidal field energy (B?) at a particu- 
lar position in space (at the base of the convection zone and at 
20°in latitude in the Northern hemisphere). Time is in years and 
the magnetic field in Gauss. Red (blue) colours indicate positive 
(negative) field values, ranging from -2 x 10^ G to 2 x 10^ G. 



The lower panel of Figure [2] shows the time evolution of the 
toroidal field energy at the base of the convection zone and at the 
latitude of 20°, inside the activity belt. We clearly note that for 
the choice of parameters we made, the solution is harmonic in 
time, no modulation of the cycle is visible. As pointed out in the 
introduction, Cha rbonneau et al. 2005l showed that a variation of 
the strength of the Babcock-Leighton poloidal source term could 
lead to a chaotic behaviour of the solutions through a sequence 
of bifurcations. In order to focus on the influence of the time 
delays on the modulation of the cycle, we chose here a value of 
Cs which lies in the range of values for which no modulation of 
the cycle can arise for the standard solution. 

4.2. Introducing moderate delays 

We now turn to investigate the effect of introducing a magnetic 
field-dependent delay in the model, modifying the surface source 
term according to Eq. |9] In particular, we focus on the modifi- 
cation of the butterfly diagram and the evolution of the toroidal 
field energy at a particular point in space. 
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layed by almost a year in this case (compared to a few hours in 
the previous calculation). We point out again that the maximum 
delay which can be reached for very weak tubes is about 6 years 
(approximately half a cycle), since a longer time for toroidal 
structures to cross the convection zone would be unrealistic. 

Figure|4]shows the butterfly diagram and the evolution of the 
magnetic energy for this case. This case is particularly interest- 
ing when we turn to consider the evolution of the toroidal energy 
with time. Indeed, again the effects on the shape of butterfly di- 
agram stay weak, even if some perturbations are clearly visible. 
In particular, close to the equator, consecutive cycles of the same 
polarity seem to connect, leading the cycle period at these lati- 
tudes to be different from the basic 11 -yr cycle still present at 
higher latitudes. We also note that the period varies from one cy- 
cle to another but this particular feature is even more visible on 
the lower panel of Figure |4] 
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Fig. 3. Same as Fig.|2]but with a delay of 14 days on 1 kG fields. 



Figure |3] shows the butterfly diagram and toroidal field en- 
ergy for a few cycles in the case where a small delay is intro- 
duced. More precisely, a delay of 14 days is imposed on 1 kG 
fields at the base of the convection zone. According to the de- 
pendence of the delay on the magnetic field intensity, it means 
that 10 kG fields will only be delayed by less than 4 hours, which 
is extremely small compared to the time scale of the magnetic 
cycle. However, even if the eff'ects are hardly visible on the but- 
terfly diagram (it is almost indistinguishable from the one of Fig. 
IS, a clear modulation of the magnetic energy appears on the 
lower panel of Fig. [3] The small modification of the butterfly di- 
agram is not surprising considering the fact that the same physi- 
cal processes still act in our model to generate and advect the 3 
components of the magnetic field. In particular, the meridional 
flow, responsible for the advection of toroidal fields at the base 
of the convection zone, still acts to produce the characteristic 
equatorward branch of sunspot migration. What is more striking 
is the effect on the cycle amplitude which becomes modulated, 
while the cycle period remains approximately constant and iden- 
tical to the standard case. To understand this particular feature, 
we will proceed to the study of a reduced model in the following 
section. But first, the effects of higher delays are analysed. 

4.3. Strong delays leading to highly perturbed solutions 

The time delay is now increased from 14 days on 1 kG fields to 
14 days on 50 kG fields, meaning that 10 kG fields will be de- 
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Fig. 4. Same as Fig. |2]but with a delay of 14 days on 50 kG 
fields. 



The major effect of the delays is again to modulate the am- 
plitude of the cycle to such a point that the weakest cycle on 
the particular period shown on Figure |4] reaches in terms of 
peak energy only 30 % of the peak energy of the strongest 
cycle. These variations are comparable to the observations of 
the monthly sunspot number available from 1750 until today 



(e.g. Hoyt & Schatten 1998) . Indeed, over this period of time, 
the weakest cycle (cycle 5, which reached its peak in 1804) has 
been shown to possess a sunspot number of about 50 whereas the 
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Strongest cycle (cycle 19, which reached its maximum value in 
1955) reached a sunspot number of more than 200. These strong 
variations are thus reproduced in this model with rather strong 
delays, but still small compared to the cycle period. 

Finally, we should stress the variety of morphologies exhib- 
ited by the different cycles. In particular, some cycles show a 
clear asymmetry between their rising and declining phases (for 
example the cycle peaking at t=895). Interestingly, we also note 
that some cycles seem to possess almost two peaks (as was ob- 
served on cycle 23) or at least show a change of slope in the 
declining phase, which makes them quite different from the stan- 
dard more "classical" cycle. 

To better understand the results of such a modulation in the 
amplitude and possibly on the period of the modelled solar cycle, 
we need to get back to a reduced model and analyse its properties 
from a dynamical systems point of view. 

5. Understanding the modulation: reduction to a 6th 
order system 

5.1. Formulation and justification of tlie reduced model 

The reduced model we decide to use does not explicitly con- 
tains time-delays as in the 2D model. It rather introduces another 
variable Q which will be delayed by time t with respect to the 
toroidal field B. Moreover, we now work in only one dimension 
in space. The set of equations we thus want to solve here is the 
following: 
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This model is significantly simplified since the physical ingre- 
dients (the rotation rate, the magnetic diffusivity, the meridional 
flow and the source term) are taken to be constants, t is still 
dependent on the inverse of the toroidal field B, leading to the 
following expression: 

To 



1 + W 



We win nevertheless show that such a simple model with a vary- 
ing time delay, even small, is able to reproduce a modulation in 
the cyclic behaviour of the solutions. 

In order to simplify even more the previous set of equations, 
we take the variables A, B and Q to be single Fourier modes in 
X, with wavenumber k, such that: 

A{x, t) - A,{t) exp(fc) 

B{x, t) - B,{t) exp(/A;jt:) 

Q{x, t) = Q,{t) expiikx) 

we then get the following set of complex ODEs: 



dt \+ iiap ' 

dB, , 

— - + ikvB, = ikOA, - nk^B, 
dt 
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dt T 
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with T = To/(l + |B/p). All the analysis of the previous set of 
ODEs will be made considering k - 1 to be the wavenumber of 
interest in our model. This system possesses 6 degrees of free- 
dom represented by the real and imaginary parts of each vari- 
ables A,, B, and Q,. 

This choice of model was motivated by the fact that we want 
to study this set of ODEs as a dynamical system. Having explicit 
delays in the equations makes it much more difficult to identify 
the sequence of bifurcations which could occur between the har- 
monic solution of Fig. |2] and the aperiodic solution of Fig. |4] 
Nevertheless, the new variable Q, introduced in our simplified 
model and which now represents the delayed toroidal field be- 
haves exactly as we want it to. To illustrate this, we solve the 
set of equations numerically for the following choice of param- 
eters: n = 100, S = 100, 77 = 10-\ V = 10-2, X = 10-4 and 
To = 2 X lO-^Bj where fii = 3 x 10^ is a normalization factor 
From now on, the values of tq will be expressed in terms of mul- 
tiples of B\. We note here that the actual time delays should not 
be directly related to the values of To, beacuse of this normaliza- 
tion factor 

These parameters are difficult to relate to the ones used in the 
full 2D model but they were chosen so that: 

- dynamo action occurs (the dynamo number is high enough), 

- the standard solution is harmonic in time (the dynamo num- 
ber is low enough so it does not lead to any modulation of 
the cycle), 

- we lie in the advection-dominated regime rather than in the 
diffusion-dominated one (the advection time 1 jkv is shorter 
than the diffusive time l/k^rj, leading to v > 77 in our case 
where A: = 1), 

- the time delay is always smaller than half a cycle period (i.e. 
we impose a maximum value for the time delay, like we did 
in the 2D model where the source was set to zero when the 
delay was too long). 
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Fig. 5. Temporal evolution of the toroidal field (plain line) and 
delayed toroidal field (dotted line) for a case where the delay is 
small enough to keep the harmonic solution stable. 



Figure|5]shows the evolution over a few cycles of the toroidal 
field %(B,) and the delayed toroidal field %(Q,). We can no- 
tice two major particularities: first, we clearly note the time shift 
(14) between the two fields, representing the time delay imposed 
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by the value of tq. In this case, the solution is harmonic and 
hence \B,\^ is constant and thus the time delay is also constant. 
Consequently, the time shift between the two fields is always 
the same along the different cycles which are represented on the 
figure. 

The second point which has to be noticed is the significantly 
lower amplitude of 'K(Q,) compared to the original toroidal field 
%{Bt). This amplitude difference is due to the fact that the evo- 
lution equation for Qt is similar to a diffusion equation. This 
property models the ohmic diffusion acting on the toroidal flux 
tubes during their rise through the convection zone. Indeed, we 
can reasonably expect the flux tubes reaching the surface to 
have a lower amplitude than at the base of the convection zone. 
Moreover, the decrease in strength will be higher when the ini- 
tial magnetic field intensity in the flux tubes is lower Indeed, 
when the flux tubes are weaker, their rise time is longer and be- 
comes comparable to the diffusive time associated to the mag- 
netic structure, which is equal to a^/rj when the rising flux tube 
has a circular section of radius a. In our reduced model, this ef- 
fect is also taken into account since the the time delay for weaker 
tubes will be higher and thus the diffusion operator will be more 
efficient in this case. 

With this choice of model and parameters, we are thus close 
to the dynamo solution calculated with the STELEM code pre- 
sented in the previous section, so that we can reasonably relate 
the results of the reduced model to the full 2D calculations. The 
major interest of this type of models is that we are now reduced 
to a 6th order system which can be analysed both analytically 
and numerically. 

5.2. Destabilization of the harmonic solution and infiuence of 
parameters 

We now proceed to the analytic study of the characteristics of 
the harmonic solution when the delay is increased. In particular, 
we investigate the influence of the model parameters on the am- 
plitude and frequency of the harmonic solution. The details of 
this analysis are presented in appendixlAl 

As expected, the dynamo number OS /rj controls the ampli- 
tude of our solution while the meridional flow speed sets the 
cycle period. This behaviour is characteristic of mean-field flux 
transport dynamo model and we are thus not going into too many 
details. The effects of an increasing delay on both the frequency 
and the amplitude of the solution is on the contrary of some in- 
terest for our particular study. 
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Fig. 6. Variation of the amplitude and the frequency of the 
toroidal field B, with the delay, for a particular case. The choice 
of parameters is the same as in the previous section, namely: 
Q = 100, S = 100, J] = 10-\ V = 10-2, A = lO-'*. 

Figure |6] shows the variation of the amplitude and frequency 
of the harmonic solution when the delay is increased when all 
the other parameters are fixed to the values quoted in the previ- 



ous section. For higher values of the delay, the solution becomes 
unstable, as we will see in the stability analysis below. The stan- 
dard "undelayed" solution is also included in these plots, for 
T = To = 0. We note that the amplitude of our solution decreases 
linearly with the delay which can be explained by the fact that 
the amplitude of 51(2,) (which represents the delayed toroidal 
field and thus is the one responsible for the regeneration of the 
poloidal field A,) is smaller than the amplitude of B,, due to the 
diffusive effects included in our equations. Hence, the amount 
of poloidal field generated by the reduced toroidal field will be 
weaker and consequently, the toroidal field produced through the 
Q-effect acting on the poloidal field will also be reduced. 

The frequency of the solution has a different behaviour when 
T is increased. Indeed, the variation is first close to linear with t, 
meaning that the cycle period is increased linearly with the delay 
but then the slope of the curve changes and the cycle is signifi- 
cantly slowed down by a variation of the delay. This feature can 
also be understood by the fact that the regeneration of poloidal 
field will take longer and longer when the delay is increased and 
thus the dynamo loop will take longer to close. 

Another interesting point here is that we can study the effect 
of the various parameters on the threshold above which the har- 
monic solution becomes unstable. We will not focus on the effect 
of the dynamo number on this threshold since it has been studied 
in similar conditions before. Instead, we fix the dynamo num- 
ber and we investigate the evolution of this threshold when the 
meridional flow speed is increased. The detailed stability anal- 
ysis is shown in appendix |B] It has to be pointed out here that 
a Floquet analysis has been applied to our equations to investi- 
gate the nature of the bifurcation undergone by the system at that 
threshold. We find that two purely imaginary conjugate Floquet 
multipliers (reflecting the behaviour of a perturbation to the peri- 
odic solution after one cycle) cross the unit circle, thus indicating 
a Hopf bifurcation. 



0.15 



0.10 



0.05 - 



0.00 



Fig. 7. Value of tq for which the harmonic solution becomes un- 
stable with respect to the normalized meridional flow speed. 



Figure [T] shows the dependence of the critical tq (i.e. the 
value for which the harmonic solution becomes unstable) with 
the meridional flow speed normalized by the magnetic diffusiv- 
ity. In our particular case where k = I (and thus where v/rj repre- 
sents the ratio between diffusive and advective times), this gives 
us some insight on the value of the delay which has to be reached 
to get a modulation, when the system tends more and more to- 
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wards an advection dominated regime. The trend of the threshold 
to decrease with increasing meridional flow speed is not surpris- 
ing if we consider we are introducing fixed values for the delays, 
not normalized to the cycle period. Indeed, as v is increased, the 
cycle period is decreased and then the absolute value for the de- 
lay to have a significant effect on the cycle is decreased. What 
is more striking is the saturation that the curve seems to reach 
when V is increased to strong values. In this case, the cycle pe- 
riod (which is roughly inversely proportional to the merid- 
ional flow speed in this type of models) is still decreased but 
the critical value for the delay seems to reach a minimum value. 
We thus conclude that the percentage of the cycle represented by 
the critical delay will increase when the meridional flow speed is 
increased. As an example, in the case where the meridional flow 
speed is quite weak (y/77 - 4), the threshold delay represents 
about 14% of the cycle period whereas at v/77 = 40, this percent- 
age reaches 30%. In this strongly advection-dominated regime, 
the meridional flow speed is probably dominating the time-scale 
of the system and as a result the effect of the delays is reduced. 

In the remainder of the paper, we focus on an intermediate 
case {v/rj = 10), where the threshold value for the delay (tq x 
3.36 X 10"^) reaches 19% of the cycle. This value can appear 
quite large to the typical values for the delay used in the full 2D 
model but it has to be pointed out that as soon as the modulation 
appears, the strongest fields reach higher values than the peak 
values of the harmonic solution. As a consequence, the delays 
for the strongest fields are shorter and reach about 8% of the 
cycle, which is closer to what was used in the 2D model. This 
remark is related to the effect of fixed versus varying delays, 
which we discuss in the conclusion. 

5.3. Further evolution 

Knowing the behaviour of the harmonic solution (which is stable 
for small values of the delay) and the influence of the parameters, 
we choose a particular set of those parameters to try to study a 
generic solution when the delay is further increased. We thus 
choose numbers which satisfy the conditions invoked in section 
15.11 and a meridional flow speed in agreement with what was 
shown in the previous section. The parameters chosen for all 
the analysis are the following: Q = 100, S - 100, 77 - 10"^, 
V = 10"^, A = 10 "^, the same as the ones used for the solution 
shown on Figure |6] 

Figure [8] shows the time series of the toroidal magnetic en- 
ergy for the solutions of the system when the delay is further 
increased to ro » 9 x 10"^. An aperiodic modulation of the solar 
cycle is here clearly visible, and the major point is that the evo- 
lution of toroidal energy is very similar to what was found in the 
full 2D system when the delay was sufficiently large (see Fig.|4|. 
Indeed, the major modulation acts on the cycle amplitude but the 
frequency also varies and periods of very low activity arise. This 
type of behaviour is very attractive when compared to the ob- 
served sunspot number which qualitatively shows the same kind 
of long-term variations. We should stress that we are still here in 
a regime where the time delays are small compared to the long- 
term modulation they seem to create. Moreover, this perturbed 
solution is obtained both through 2D mean-field simulations and 
through a reduced set of ODEs solving the Babcock-Leighton 
dynamo problem. These results are thus of particular interest as 
a new origin of modulation seems to appear: the time delays due 
to the buoyant rise of flux tubes from the base of the convec- 
tion zone to the solar surface to create active regions. We now 
need to understand what happens to the system between the first 
bifurcation and this aperiodic more interesting solution. 
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Fig. 8. Time evolution of the toroidal magnetic energy with a 
value of the delay of tq = 9 x 10"^. Aperiodic modulation of the 
cyclic pattern is now fully developed. 



We thus turn to investigate in more details which sequence 
of bifurcations lead to such complicated behaviour. 

6. Beyond the first bifurcation: reduction to a 5th 
order system 

To gain some insight on the behaviour of the system after the 
first bifurcation, we show that we can further reduce our prob- 
lem from a 6th order to a 5th order system. Indeed, it is rather 
impractical to study the precise behaviour of the system after the 
first bifurcation. Reducing our system even more will help us 
achieve this goal. 

6.1. Formulation of the model 

Noticing that our system possesses a phase-invariance, i.e. a fur- 
ther symmetry, we follow the procedure of Weiss et al. (1984) 
and |Jones et al. (1985)] and make a change of variables which 
will remove this symmetry. 

We express our variables A,, B, and Q, using polar coordi- 
nates and introduce four new variables as shown 
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where y and z are complex numbers and and p are real. The 
real numbers p and 9 here represent respectively the modulus and 
argument of the basic toroidal field B,. We can now reintroduce 
these expressions in equations [T3l fT4l and [TSJ to get, after some 
algebra, the new set of ODEs 
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We then get a 5th order system since y and z are complex 
but p is real. We find that 6, which is the phase of the "unde- 
layed" toroidal field, does not appear explicitly in this new set of 
equations but can be obtained from 



e = np%(y)-vp 



(19) 



The interesting outcome of reducing our system from 6 de- 
grees of freedom to 5 is that the limit cycle of the higher-order 
system now corresponds to a stationary solution in our reduced 
model and that a torus now corresponds to a limit cycle. 

Solving this new system numerically will now enable us to 
analyze the further evolution of our solution when the delay is 
increased. In particular, we show that we can now easily exhibit 
the sequence of bifurcations leading to the aperiodic modulation 
of Fig. El 

6.2. From the harmonic to the aperiodic solution: a sequence 
of bifurcations 

We follow the evolution of our system from the first bifurcation 
at To - 3.36 X 10"^ to the aperiodic solution at tq = 9 x 10"^. 
First, we can relate again these values to some physically mean- 
ingful time-delays, keeping in mind that the system has been 
simplified and that a particular choice of parameters was made. 
For the stronger fields, a value of tq ~ 3.36 x 10"^ would then 
correspond to a delay of approximately a tenth of the cycle du- 
ration and To ~ 9 X 10"^ a third of the cycle period. These delays 
may look higher than what was calculated in the 2D model but 
we showed that a modification of the parameters (and especially 
of the meridional flow speed) could move the threshold above 
which the modulation starts to appear. 

The Floquet multipliers for this system are computed in or- 
der to study the stability of our periodic solution. We find that 
the solution of our new reduced model stays on a periodic orbit 
for a large range of values for the delay, i.e. from 3.36 x 10"^ to 
about 5.43 x 10"^. We indeed find that at this value for the param- 
eter To, one of the real Floquet multipliers crosses the unit circle 
through the -1 point, indicating a period doubling bifurcation. 
This new 2-periodic solution then persists for quite a large range 
of the control parameter until we reach a new period doubling 
bifurcation. 

Figure |9] shows the evolution of the system from the peri- 
odic solution resulting from the first Hopf bifurcation at tq ~ 
3.36 X 10"^ until we reach the aperiodic modulations of the cy- 
cle. We present on this figure a bifurcation diagram computed 
by determining the maximum values of p for each value of 
To > 4 X 10"^ (when the first period doubling has already oc- 
cured). The initial branch existing at the beginning of the bifur- 
cation diagram thus indicates a periodic orbit. Moreover, phase 
portraits at key values of the parameter are shown. They were 
chosen to stress the sequence of period doubling bifurcations ex- 
hibited by the set of ODEs. Panel b) shows the modification of 
the phase portrait of panel a) after the first period doubling has 
occured. We are then in the presence of a 2-periodic solution: ev- 
ery 2 cycles, the value of p reaches the same maximum value and 
this behaviour persists for as long as the model has been calcu- 
lated, i.e. for thousands of cycles. At To ~ 7.15x10"^, the system 
goes through another bifurcation. The period is doubled again, 
which is visible both on the bifurcation diagram (where the num- 
ber of points is doubled at this value) and on panel c). Indeed, the 
curve on the phase portrait closes on itself after 4 cycles now in- 
stead of 2 for the previous values of to. Period doublings happen 
again at to ^ 7.8 x 10"^ and To ~ 7.88 x 10"^ where the so- 
lution becomes first 8-periodic and then 16-periodic. The phase 



portrait on panel d) shows the time needed to close the curve is 
doubled at to ~ 7.8 x 10"^ compared to the solution of panel c). 
But obviously, since the curve is still closed, the system is still 
periodic. 

It then becomes difficult to identify further period doublings 
but they may occur in very narrow windows for the parameter. 
The main point is that the system seems to reach a very com- 
plicated behaviour when the value of to approaches 8.5 x 10"^. 
Panel e) shows the phase portrait for to - 8.8 x 10"^. The phase 
space is almost entirely covered by the solution and any value 
can be reached by the peaks in the temporal evolution of p, as 
shown on the bifurcation diagram. At this stage, a clear periodic- 
ity is difficult to identify and we recover the strongly modulated 
solution shown on Fig . [8]for a value of the parameter of 9 X 1 0~^ . 

Consequently, our system seems to undergo a sequence 
of period doublings leading to a chaotic behaviour. The so- 
lution may well be modified for higher values of the time- 
delays but which would not be realistic anymore for the 
physical process we are modelling. Similar cascades of pe- 
riod doublings have been found in the Lorenz equations 
([S parrow 1982[ ) and in other systems of ordinary and partial dif- 
ferential equations, including those governing the solar dynamo 
(.Moore et al.T983l IKnobloch & Weiss 19831 IWeiss et al. 19841 
IKnobloch et al. 1998I ). However, the nonlinearities at the origin 
of the sequence of bifurcations in those models were always 
related to the dynamical feedback of the Lorentz force on the 
flow. The striking result in this work is that modulation of the 
cycle amplitude can also arise in models where the only nonlin- 
earities are the quenching term in the Babcock-Leighton source 
and the magnetic field-dependent time delays. It has here to be 
noted though that the time delays also appear in the quenching 
term (see Eq.|9]l, possibly amplifying the nonlinear effects in the 
model. 

We conclude that time delays due to flux tubes rising from 
the base of the convection zone to the surface seem to create 
a strongly modulated cycle in Babcock-Leighton dynamo mod- 
els. A sequence of period doublings leading to chaos has clearly 
been identified here thanks to our 5th-order system. Indeed, even 
though it is also obviously present in the 6th order system, it 
is much more difficult to exhibit. The magnetic cycle result- 
ing from these models which take into account the rise time of 
toroidal fields thus seems to be closer to the actual behaviour of 
the solar cycle. 

7. Discussion and conclusion 

In this work, we have introduced a more physically accurate 
model for the poloidal field source term in Babcock-Leighton 
flux transport dynamo models. To do so, some results from 3D 
MHD calculations of rising flux tubes from the base of the con- 
vection zone to the surface were reintroduced in 2D mean-field 
models. In particular, the rise time of flux ropes, which was as- 
sumed to be negligible in previous simulations, was here taken 
into account by introducing magnetic energy-dependent time de- 
lays in the BL source term. Through the use of a combination 
of analytical and numerical techniques applied to an adapted 
reduced system of non linear equations, we show that the sys- 
tem exhibits a sequence of bifurcations leading to a chaotic be- 
haviour. The time evolution of the solution at that stage is qual- 
itatively very similar to the full 2D model and more importantly 
to the actual solar activity, with strong modulation especially on 
the cycle amplitude. 

Assuming the rise of flux tubes from the base of the con- 
vection zone and the surface to be instantaneous and thus the 
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Fig. 9. Evolution of the system after tq = 4 x 10 ^: bifurcation diagram and phase portraits in the (p, 5^(pz)) plane for values of the 
control parameter tq close to those for which successive period doublings occur. 



time delays due to magnetic buoyancy to be negligible can sound 
reasonable since these delays are weak compared to the other 
time-scales of the system. In particular, we can expect the time 
delay due to the advection by the meridional flow (which has a 
characteristic time-scale of a few solar cycles, see for example 
Charbonneau & Dikpati 2000) l to be much more eff'ective on the 
behaviour of the dynamo-generated magnetic field. As a con- 
sequence, it is mainly the eff'ect of meridional flow which has 
been tested in previous models (e.g. Wil mot-Smith et al. 2005I) . 
However, we showed here that the time delays due to flux tubes 
rising have in fact a significant efifect to modulate the cycle on 
time scales very large compared to the delays. The main reason 
for this behaviour is that the delays are not fixed but depend on 
the modulated toroidal energy itself. To illustrate this explana- 
tion, we show on Figure [10] the time evolution of toroidal fields 
(delayed and undelayed) for a case similar to what we showed in 
the core ofthe paper (a |Bp-dependentdelay with To - 3.4x10"^, 
i.e. just after the Hopf bifurcation) and the same evolution when 



the delay is fixed to the mean value reached by the delay in the 
previous run. 



On this figure, both the basic toroidal field and the delayed 
field are shown, stressing the time shift existing between the two. 
The striking result is that a long-term modulation appears only 
in the case of a varying delay. In the fixed-delay case, although 
the time-shift between the basic toroidal field and its delayed 
counterpart is quite significant compared to the cycle period, no 
modulation is created. On the contrary, in the varying-delay case, 
even though the strongest fields are almost not affected by the 
time-shift, a strong modulation on the amplitude is built up. As 
a result, the fact that small delays can aff'ect the long-term evo- 
lution of the dynamo cycle seems to be linked to the variability 
of these delays and more specifically to their dependence on the 
magnetic field strength. Since magnetic buoyancy is much more 
efficient in the solar interior for strong toroidal field structures, 
these kind of delays are likely to appear in reality. Consequently, 
a new origin of the solar cycle modulation may have been re- 
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Fig. 10. Time evolution of the toroidal fields amplitude (unde- 
layed in plain line and delayed in dotted line) for a delay depen- 
dent on |Bp just after the Hopf bifurcation (upper panel) and for 
a delay fixed to the mean value of the previous run (lower panel). 



vealed here, thanks to the input of 3D MHD models into 2D 
mean-field dynamo calculations. 

Only a part of the results of 3D calculations were reintro- 
duced here but other effects can also be taken into account. In 
particular, hoop stresses and the Coriolis force are known since 
the first thin flux tubes simulations to deflect the trajectory of 
buoyantly rising magnetic fields inside the convection zone, this 
effect being stronger when the fields are less intense. This feature 
could also be introduced in BL models by modifying the source 
term for poloidal field. In particular, the source term is likely to 
get a contribution from the same latitude a the initial position of 
the flux tubes when those are strong enough and from a lower 
latitude when the flux tubes have been sufficiently influenced by 
the Coriolis force and the hoop stresses to rise parallel to the ro- 
tation axis. This may result in a different shape for the butterfly 
diagram with slightly more activity at higher latitudes, provided 
that weak magnetic structures are assumed to be able to reach 
the surface and create active regions. This could lead to recon- 
sider the correspondence which is often made between toroidal 
magnetic field generated in 2D dynamo models at the base of 
the convection zone and the actual sunspot migration observed 
at the solar surface. 

The variations of the solar cycle period are also known to be 
significant. For instance, cycle 23 has been significantly longer 
than the previous ones, with a duration of about 13 years (1996- 
2009). In the models computed in this work, the modulation of 



the solar activity is particularly visible on the cycle amplitude 
and is much less obvious on the period. However, the influence 
of rising flux tubes on the meridional flow was not taken into ac- 
count here. Since the cycle period of the solutions resulting from 
these flux-transport models is very sensitive to the amplitude of 
the meridional circulation, a modification of the flow by rising 
flux tubes is very likely to modify the frequency of the dynamo 
solution and thus to also introduce a modulation on the cycle 
period. This feature is currently being investigated. 

A considerable step forward would obviously be to develop 
a self-consistent global model with buoyant toroidal structures 
built up and making their way from the base of the convection 
zone where they become unstable to the photosphere where they 
create active regions. Unfortunately, this has not been achieved 
yet due to numerous physical and numerical difficulties. In the 
mean time, we have shown here that inputs from 3D MHD mod- 
els simulating a particular step of the dynamo cycle can sig- 
nificantly improve 2D mean-field calculations. In particular, we 
have shown that they can even help reproducing a feature of the 
solar cycle (its variability) which was not present in the standard 
model before. We will continue to develop these ideas in future 
work. 
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Appendix A: Amplitude and frequency of the 
harmonic solution 

We can further expand our harmonic solution on a single Fourier 
mode in time as: 

A,{t) - Ao exp(;wf) 
B,{f) - Bq exp(/wf) 
Qtif) = Go exp(/wO 

Reintroducing these expansions in equations [13] [14] and [15] 
will give us a way to calculate the frequency u and the amplitude 
of the harmonic solution |B()| or equivalently \Qq\ with respect to 
all the parameters. The complex dispersion relation we get with 
^ = 1 is the following 



{iu + iv + rj) {i(jjT+l) 



ins 



1 + ^ieo 



(A.l) 



Taking the real and imaginary parts of this relation dispersion 
give us the following equations 



w^(l +2t]t) + 2cov{1+t?]) + v^ -T]'^ ^0 



ll) T + 2aj VT + cl>(v t — rj t — Irf) 



QS 



1 + ^ieoi' 



(A.2) 



+ 277v(A.3) 



where 



T = To/(l + IBoP) 



where ai, a^, j6i, fii, y\ and 72 represent the coefficients of the 
perturbed fields, p represents the complex growth rate of the per- 
turbation and the symbol • stands for the complex conjugate. 

Substituting these expressions in equations [T3] [14] and [15] 
give us a system of 6 equations relating the coefficients of the 
perturbations, the growth rate and all the parameters of the sys- 
tem. 
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|BoP and |2oP being related by the following expression 

\Bq\H\ + |Bol¥ 



■ySi 



Igor = 



co^Tf^ + (l+\Bo\^) 



2\2 



A non-trivial solution for this system exists only if the de- 
terminant is zero, we thus have to calculate the determinant and 
find the roots of the 6th order polynomial, in terms of the growth 
rate p. The harmonic solution loses stability at a Hopf bifurca- 
tion when the real part of one the roots becomes positive. 



Appendix B: Stability of the harmonic solution 

We can focus on the stability of this periodic solution by per- 
turbing it and studying the growth rates of the perturbations. To 
do so, we perturb the solution as follows 



A = A, (l-HQ'ie'"-H «*£''*') 
B^B,(l+/3ie'" +/3*eP*') 
Q^ Qr (I + 71 eP'+y*eP*') 



(B.l) 
(B.2) 
(B.3) 



